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Abstract. A wide variety of articles, starting with the famous paper [TT], is devoted 
to the uniqueness question for the semilinear elliptic boundary value problem — Ati = 
Xu + u p in il, u > in ft, u — on dfl, where A ranges between and the first Dirichlet 
Laplacian eigenvalue. So far, this question was settled in the case of f2 being a ball 
and, for more general domains, in the case A = 0. In |16j . we proposed a computer- 
^vjj | assisted approach to this uniqueness question, which indeed provided a proof in the case 

il = (0, l) 2 , and p = 2. Due to the high numerical complexity, we were not able in (16) 
to treat higher values of p. Here, by a significant reduction of the complexity, we will 
^1 . prove uniqueness for the case p — 3. 
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Dedicated to the memory of Wolfgang Walter 

1. Introduction 

The semilinear elliptic boundary value problem 



- A« = /(«) in fi, u = 0onffl (1.1) 

has attracted a lot of attention since the 19th century. Questions of existence and mul- 
tiplicity have been (are still being) extensively studied by means of variational methods, 
fixed-point methods, sub- and supersolutions, index and degree theory, and more. 
In this article, we will address the question of uniqueness of solutions for the more special 
(T) I problem 

OV / 

I -Am = Xu + u p in f2 

I u > in Q (1.2) 

© \ u = on dtt 

t— i ; v 

where A ranges between and Ai(f2), the first eigenvalue of the Dirichlet Laplacian. It 
has been shown in a series of papers [19], [29], [28], [T], [2] that the solution of ( II. 2p is 
indeed unique when Q is a ball, or when Q is more general but A = ([30], [12], [9], [TO]). 
We will concentrate on the case where Q = (0, l) 2 and p = 3, and prove that uniqueness 
holds for the full range [0, Ai(fi)) of A. Thus, our paper constitutes the first uniqueness 
result for this situation. More precisely we prove 

Theorem 1.1. Let Q be the unit square in M?, Q = (0, l) 2 . Then the problem 

—Am = Xu + m 3 in Q 

u > in Q (1.3) 

u = on dQ 

admits only one solution for any X G [0, Ai(fi)). 
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Remark 1.1. a) A simple scaling argument shows that our uniqueness result carries over 
to all squares Qi := (0,Z) 2 (and thus, to all squares in M 2 ): If u is a positive solution of 
—Am = Xu + u 3 in Qi, u = on dQ[, for some A G [0, Ai(f2;)), then v(x, y) := lu(lx, ly) is 
a solution of (TO} for A = AZ 2 G [0, Ai(fi)). 

b) Since we also show that the unique solution in the square is nondegenerate, by a result 
of [I0] we deduce that the solution is unique also in domains "close to" a square. 

c) Finally we observe that having shown in [16] (case p = 2) and in this paper (case 
p = 3) that the unique solution is nondegenerate then uniqueness follows also for other 
nonlinearities of the type Xu + u p for p close to 2 and 3. Indeed, by standard arguments 
(see for example [9]) nonuniqueness of positive solutions in correspondence to sequences 
of exponents converging to 3 (resp. to 2) would imply degeneracy of the solution for p = 3 
(resp. p = 2). 

Our proof heavily relies on computer- assistance. Such computer-assisted proofs are re- 
ceiving an increasing attention in the recent years since such methods provided results 
which apparently could not be obtained by purely analytical means (see [6], [5], [21], [T7] . 

M)- 

We compute a branch of approximate solutions and prove existence of a true solution 
branch close to it, using fixed point techniques. By eigenvalue enclosure methods, and 
an additional analytical argument for A close to Ai(f2) we deduce the non-degeneracy of 
all solutions along this branch, whence uniqueness follows from the known bifurcation 
structure of the problem. 

In [16] we give a general description of these computer-assisted means and use them to 
obtain the desired uniqueness result for the case Q = (0, l) 2 ,p = 2. To make the present 
paper dealing with the case p = 3 more self-contained, we recall parts of the content of 
[16] here. We remark that the numerical tools used in [T6] turned out not to be sufficient 
to treat the case p = 3. Now, by some new trick to reduce the numerical complexity, we 
are able to handle this case. 



2. Preliminaries 

In the following, let Q = (0, l) 2 . We remark that the results of this section can be carried 
over to the more general case of a "doubly symmetric" domain; see [16] for details. 
First, note that problem (11. 2p can equivalently be reformulated as finding a non-trivial 
solution of 

-Au = Xu + \u\ p in Q , > 

u = o on on, { ' 

since, for A < Ai(f2), by the strong maximum principle (for —A — A) every non-trivial 
solution of (12. ip is positive in Q. In fact, this formulation is better suited for our computer- 
assisted approach than (11. 2p . 

As a consequence of the classical bifurcation theorem of j2S] and of the results of [H] the 
following result was obtained in [20] : 

Theorem 2.1. All solutions u\ of (ll.2p lie on a simple continuous curve T in [0, Ai(f2)) x 
C 1,a (Q) joining (Ai(O),0) with (0,Uq), where Uq is the unique solution of (II. 2p for A = 0. 

We recall that the uniqueness of the solution of (11.21) for A = was proved in [TU] and [S]. 
As a consequence of the previous theorem we have 
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Corollary 2.1. If all solutions on the curve T are nondegenerate then problem (\1.2\i 
admits only one solution for every A G [0, Ai(fi)). 

Proof. The nondegeneracy of the solutions implies, by the Implicit Function Theorem, 
that neither turning points nor secondary bifurcations can exist along T. Then, for every 
A G [0, Ai(fi)) there exists only one solution of (11.21) on T. By Theorem 12.11 all solutions 
are on T, hence uniqueness follows. □ 

Theorem 12.11 and Corollary 12.11 indicate that to prove the uniqueness of the solution of 
problem (11.21) for every A G [0, Ai(fi)) it is enough to construct a branch of nondegenerate 
solutions which connects (0,Wo) to (Ai(f2),0). This is what we will do numerically in the 
next sections with a rigorous computer-assisted proof. 

However, establishing the nondegeneracy of solutions u\ for A close to Ai(fi) numerically 
can be difficult, due to the fact that the only solution at A = Ai(f2), which is the identi- 
cally zero solution, is obviously degenerate because its linearized operator is L = — A — Ai 
which has the first eigenvalue equal to zero. The next proposition shows that there exists 
a computable number A(f2) G (0, Ai(f2)) such that for any A G [A(O) , Ai(f2)) problem (11.21) 
has only one solution which is also nondegenerate. Of course, from the well-known results 
of Crandall and Rabinowitz, [7J [8], one can establish that for A "close to" Ai, all solutions 
u\ are nondegenerate. However, in order to complete our program, we need to calculate 
a precise and explicit estimate of how close they need to be. This allows us to carry out 
the numerical computation only in the interval [0, A(fi)] as we will do later. 

Let us denote by Ai = Ai(fi) and A 2 = A 2 (fi) the first and second eigenvalue of the 
operator —A in Q with homogeneous Dirichlet boundary conditions. We have 

Proposition 2.1. // there exists A G (0, Ai) and a solution u\ of (II. 2p with A = A such 
that 

then 



P J VAi 



MU<('^^T\ (2.3) 



V 

and U\ is non- degenerate, for all solutions U\ of (II. 2p belonging to the branch r 2 C T 
which connects (\,ux) to (Ai,0). 

(Recall that V is the unique continuous branch of solutions given by Theorem 12.11 ) 
Proof, see [16] 

Corollary 2.2. If on the branch V there exists a solution u^, A G (0, Ai) such that: 

i) on the sub-branch Ti connecting (0, Uq) with (A, u^) all solutions are nondegenerate 

and 

then all solutions of (\l.2\i are nondegenerate, for all A G (0, Ai) ; and therefore problem 
(II. 2p admits only one solution for every A G [0, Ai(f2)). 
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Proof. We set r = I\ U T 2 with Ti connecting (0,w ) to (\,Ux)- On Ti we have that 
all solutions are nondegenerate by i). On the other hand the hypothesis ii) allows us to 
apply Proposition 12.11 which shows nondegeneracy of all solutions on IY Hence there is 
nondegeneracy all along T so the assertion follows from Corollary 12.11 □ 

The last corollary suggests the method of proving the uniqueness through computer as- 
sistance: first we construct a branch of nondegenerate "true" solutions near approximate 
ones in a certain interval [0, A] and then verify ii) for the solution u\. Note that the esti- 
mate f )2.4p depends only on p and on the eigenvalues Ai and A 2 of the operator —A in the 
domain Q. So the constant on the right-hand side is easily computable. When Q is the 
unit square which is the case analyzed in the next sections, the estimate (12 .4p becomes: 

37i 2 \^ / A \ ^ /3A\^ 



because Ai = 2tt and A2 = 5% . 



Fixing p = 3 we finally get the condition 



Klloo<y^. (2.5) 



3. The basic existence and enclosure theorem 

We start the computer-assisted part of our proof with a basic theorem on existence, local 
uniqueness, and non- degeneracy of solutions to problem ( 12. ip . assuming p = 3 now for 
simplicity of presentation. In this section, the parameter A G [0, Ai(f2)) is fixed. 
Let Hq (f2) be endowed with the inner product (u, v) H i := (Vm, Vt> )l 2 +cr(u, v) 12; actually 
we choose a = 1 in this paper, but different (usually positive) choices of a are advantageous 
or even mandatory in other applications, whence we keep a as a parameter in the following. 
Let H~ l {VL) denote the (topological) dual of Hq(Q), endowed with the usual operator 
sup —norm. 

Suppose that an approximate solution u\ G Hi (fi) of problem (12. ip has been computed 
by numerical means, and that a bound 5\ > for its defect is known, i.e. 

\\ - Au x - \uj x - \u\\ 3 \\h-^ < $\, (3.1) 

as well as a constant Kx such that 

\\v\\ Hh < Kx\\L { x^)[v}\\h^ for all v G H^Q). (3.2) 

Here, £(a,u! A ) denotes the operator linearizing problem (12.11) at ux] more generally, for 
(A,w) 6lx Hq(Q), let the linear operator L( X , U ) '■ Hl(Q) — > i/ _1 (fi) be defined by 

L { x,u)[v] := -Av -Xv- 3\u\uv (v G F x (O)). (3.3) 

The practical computation of bounds d~x and Kx will be addressed in Sections 6, 7 and 8. 
Let C4 denote a norm bound (embedding constant) for the embedding Hl(Q) L 4 (f2), 
which is bounded since Q C IR 2 . C4 can be calculated e.g. according to the explicit 
formula given in [231 Lemma 2]. Finally, let 

7 := 3C 4 3 . 
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In our example case where Q = (0, l) 2 , the above-mentioned explicit formula gives (with 
the choice a := 1) 

3V2 ( 1 



7 4(vr2 + l) 3 / 4 V 5 
Theorem 3.1. Suppose that some a x > exists such that 

6x<^-i a l(\\Lo x \\ L . + C 4 a x ) (3.4) 

and 

2K xl a x (\\u x \\ L 4 + C 4 a x )<l. (3.5) 

Then, the following statements hold true: 

a) (existence) There exists a solution u x G Hq(Q) of problem (12. ip such that 

\\u x - uj x \\ h i < a x . (3.6) 

b) (local uniqueness) Let r\ > be chosen such that (13. 5p holds with a x + r\ instead 
of a x . Then, 



u e HliVt) solution of |HD 

\\u-uj x \\ h i <a x + rj 

c) (nondegeneracy) 

u e Hl(tt) 
\\u - u x \\ H i < a x 

whence in particular £(a,« a ) is bijective (by (13. 6p ). 

For a proof, see 



u = u x . (3.7) 



l =>■ L( X)U ) : Hq(£1) H 1 (Vt) is bijective, 



(3.8) 



Corollary 3.1. Suppose that (13. 4p and (13. 5p no/d ; and m addition that \\u x \\ H i > a x . 
Then, the solution u x given by Theorem 3.1 is non-trivial (and hence positive). 

Remark 3.1. a) The function ijj{a) := jr- — lo^iWus-^n + C 4 a) has obviously a positive 

maximum at a = ^=r ^^/II^aII^ + — \\^x\\l^j > an d the crucial condition (13 .4p requires 
that 

4C 4 + 7^A||wA|li 4 



5\ < 



K x {^ 1 K x { 1 K x \\u x f LA + 3C 4 ) + 1 K x \\u x \\ Li ) 

1 



(V^a^IMI^ + 3C 4 ) + >yK x \\u x \\ L A + 6C 4 ) ' 

(3.9) 

i.e. 5a has to be sufficiently small. According to (13. ip . this means that u x must be 
computed with sufficient accuracy, which leaves the "hard work" to the computer! 
Furthermore, a "small" defect bound 5 X allows (via (13 .4p ) a "small" error bound a x , if 
K x is not too large. 

b) If moreover we choose the minimal a x satisfying (13. 4p . then the additional condition 
(13. 5p follows automatically, which can be seen as follows: the minimal choice of a x shows 



that a\< a. We have 
2Kx^a(\\u\\\ L A + C 4 a) = 

1 o ; ° 4 <1 (3-10) 

3C 4 + 2 1 K x \\tu x \\ 2 Li + 2^7^(7^11^11^ + 3C 4 )||o;a||l4 

and thus condition (I3.5P is satisfied. 

Since we will anyway try to find ot\ (satisfying (13. 4p ) close to the minimal one, condition 
(13. 5 p is "practically" always satisfied if (13 .4p holds. (Nevertheless, it must of course be 
checked.) 

4. The branch (u\) 

Fixing some A G (0, Ai(O)) (the actual choice of which is made on the basis of Proposition 
12. 1( see also Section 5), we assume now that for every A G [0, A] an approximate solution 
u\ G Hq(Q) is at hand, as well as a defect bound 5\ satisfying (13. ip . and a bound K\ 
satisfying (13. 2p . Furthermore, we assume now that, for every A G [0, A], some a\ > 
satisfies (13 ,4p and (I3.5p . and the additional non-triviality condition > a A (see 

Corollary 3.1). We suppose that some uniform (A-independent) rj > can be chosen such 
that (13. 5p holds with ci\ + 77 instead of ci\ (compare Theorem 13.11 b)). Hence Theorem 
13.11 gives a positive solution u\ G Hq{VL) of problem (12. ip with the properties (13. 6p . (13.71) . 
(E31), for every A G [0, A]. 

Finally, we assume that the approximate solution branch ([0, A] — > Hq (Q), A 1— > U\) is 
continuous, and that ([0, A] — > M, A (->■ a\) is lower semi-continuous. 

In Sections 6, 7 and 8, we will address the actual computation of such branches (uj\), (#a), (-^a) 
(«a)- 

So far we know nothing about continuity or smoothness of ([0,A] — > Hq(Q), A i-> u\), 
which however we will need to conclude that (u\) Xe [o,\] coincides with the sub-branch Ti 
introduced in Corollary 12.21 

Theorem 4.1. The solution branch 

[0,A] H^(Q) 
A i—)- «a 

is continuously differentiable. 
Proof: The mapping 

T f r x H%(n) -)> iy-^o) 

\ (A, tt) h-> — Am — Am — |w| 3 

is continuously differentiable, with (dJ-/du)(X,u) = £(a,«) (see (13.31) ). and J-(\,U\) = 
for all A G [0, A]. Using the Mean Value Theorem one can show that I/(a,«) depends indeed 
continuously on (A,w); see [TB], Lemma 3.1] for details. 

It suffices to prove the asserted smoothness locally. Thus, fix Ao G [0, A]. Since £(a ,m a ) i s 
bijective by Theorem 13.11 c). the Implicit Function Theorem gives a C 1 -smooth solution 
branch 

(X -e,X + e) -)• H*(n) 
A 1 — y zt\ 

to problem (12. ip . with ii\ = u\ . By ( 13. 6p . 

II^ao - u Xo \\m < «A • (4.1) 
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Since u\ and u\ depend continuously on A, and a\ lower semi-continuously, (14. ip implies 

\\u\ - waIIjj-i < «a + V (AG [0, A] n (A - e, A + e)) 

for some e G (0, e). Hence Theorem 13. li b) provides 

u x = u x (A G [0, A] fi (A — e, A + §)) , 

implying the desired smoothness in some neighborhood of Ao (which of course is one-sided 
if A = or A = A). □ 

As a consequence of Theorem I4.1[ (wa)a6[o,a] is a continuous solution curve connecting 
the point (0, tt ) with (A, Ux), and thus must coincide with the sub-branch I\, connecting 
these two points, of the unique simple continuous curve T given by Theorem 12.11 Using 
Theorem 13.11 c). we obtain 

Corollary 4.1. On the sub-branch F\ ofT which connects (0, uq) with (A, u^), all solutions 
are nondegenerate. 

Thus, if we can choose A such that condition (12. 4p holds true, Corollary 12.21 will give the 
desired uniqueness result. 

5. Choice of A 

We have to choose A such that condition (12 .4p is satisfied. For this purpose, we use 
computer-assistance again. With xm denoting the intersection of the symmetry axes of the 
(doubly symmetric) domain f2, i.e. xm = (§, |) for = (0, l) 2 , we choose A G (0, Ai(fi)), 
not too close to Ai(fi), such that our approximate solution u\ satisfies 

with "not too small" difference between right- and left-hand side. Such a A can be found 

within a few numerical trials. 

Here, we impose the additional requirement 

cox £ H 2 {Q) fl Hq{VI) , (5.2) 

which is in fact a condition on the numerical method used to compute U\. (Actually, 
condition (15.21) could be avoided if we were willing to accept additional technical effort.) 
Moreover, exceeding (13. ip - we will now need an L 2 -bound 8\ for the defect: 

||-Ac^-Ac^-|u; x | 3 || L2 <^. (5.3) 

Finally, we note that Q is convex, and hence in particular if 2 -regular, whence every 
solution u G i?o(0) of problem (JZTj) is in if 2 (fi). 

Using the method described in Section 3, we obtain, by Theorem 13. II a), a positive solution 
U\ G H 2 (tt) n #o(fi) of problem (JUj) satisfying 

\W\-Ux\\hi < «a , (5-4) 
provided that (13. 4p and (13.51) hold, and that || jy 1 > ct\. 

Now we make use of the explicit version of the Sobolev embedding H 2 (Q) C(U) given 
in [21]. There, explicit constants C , Ci, C 2 are computed such that 

IMloo < C'oII^IIl 2 + C'i||Vu|| i 2 + C^HuzrrlU 2 for all u G H 2 (fl) , 
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with ||m xx ||l2 denoting the L 2 -Frobenius norm of the Hessian matrix u xx . E.g. for Q = 
(0, l) 2 , [2TJ gives 

C = 1, C x = 1-1548 • J- < 0.9429, C 2 = 0.22361 • J— < 0.1764 . 

V 3 V 45 

Moreover, ||u xx ||l2 < ||Ait||x,2 for u G H 2 (Q) n Hq(Q) since Q is convex (see e.g. [14"]). 
Consequently, 

||«a - w aIU < C \\ux - uj-\\\ l z + Ci\\u- X - u x \\ H i + C 2 \\Au x - Ac^|| L 2. (5.5) 
To bound the last term on the right-hand side, we first note that 



N 3 - KI 3 IU 2 



K + t( u x - Wa)|(wa + t(«A - u- x ))dt ■ (u x - U) X ) 



< 3 / \\\u x + t(u x - u x )\ ■ \u x - u x \\\ L 2dt 



L 2 



< 3 / + -Wa)|| L 8||ua -w x || L odt 

< 3 / (HwaIIlb + tC 6 a A ) 2 ^ • Cg^a 



3C 6 M|wa||£b + C 6 ||a;A||i6Q;A + ^ C e a ~\j 



(5.6) 



(5.7) 



using (15 .4p and an embedding constant Cg for the embedding Hq(Q) > L 6 (Q) in the last 
but one line; see e.g. [231 Lemma 2] for its computation. Moreover, by (12.11) and (I5.3p . 

\\Au x - Au x \\ L 2 < 6- x + \\\ u - x -uj- x \\ L2 + \\\u x \ 3 -\u x \ 3 \\ L 2. (5.8) 

Using (15.41) - (15.81) . and the Poincare inequality 

)u\\ L 2 < \\u\\ H i (u e H*(ty), (5.9) 



a 



we finally obtain 



m-wx\ 



< 



Co + XC 2 



! (~\ s ^C'.iC'o ( ||u: A || 2 6 + C 6 ||u; A || L 6a A + \ C % a l 



■a x + C 2 8 x , (5.10) 



and the right-hand side is "small" if a x and S x are "small" , which can (again) be achieved 
by sufficiently accurate numerical computations. 

Finally, since 

U X (X M ) < U\(xm) + \\u x - WaIIoo, 

(15.1 0p yields an upper bound for u x {xm) which is "not too much" larger than u x (xm)- 
Hence, since u x (xm) = H^aIIoo by [TTj . condition (12. 4p can easily be checked, and (15. ip 
(with "not too small" difference between right- and left-hand side) implies a good chance 
that this check will be successful; otherwise, A has to be chosen a bit larger. 
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6. Computation of uj\ and 5\ for fixed A 

In this section we report on the computation of an approximate solution uj\ G H 2 (Q) D 
Hq(Q) to problem ( 12. ip . and of bounds 5\ and K\ satisfying (13. ip and ( 13. 2p . where 
A G [0, Ai(f2)) is fixed (or one of finitely many values). We will again restrict ourselves to 
the unit square Q = (0, l) 2 . 

An approximation u\ is computed by a Newton iteration applied to problem ( 12. ip . where 
the linear boundary value problems 

L {x ^ )) [v n ] = Au^ + Xut ) + \ujt^ (6.1) 
occurring in the single iteration steps are solved approximately by an ansatz 

N 

v n (xi,X2) = aff sm(i7TXi) sin( j-rrx 2) (6.2) 

i,i=l 

and a Ritz-Galerkin method (with the basis functions in (I6.2p ) applied to problem (16.11) . 
The update Co^ n+1 ^ := u^' + v n concludes the iteration step. 

The Newton iteration is terminated when the coefficients affi in (16. 2p are "small enough", 
i.e. their modulus is below some pre-assigned tolerance. 

To start the Newton iteration, i.e. to find an appropriate o;^ '' of the form (16. 2p . we 
first consider some A close to Ai(fl), and choose u^\xi,X2) = a sin(7ra;i) sin(7rx 2 ); with 
an appropriate choice of a > (to be determined in a few numerical trials), the Newton 
iteration will "converge" to a non-trivial approximation wW, Then, starting at this value, 
we diminish A in small steps until we arrive at A = 0, while in each of these steps the 
approximation computed in the previous step is taken as a start of the Newton 
iteration. In this way, we find approximations U\ to problem (12.11) for "many" values of 
A. Note that all approximations u\ obtained in this way are of the form (16. 2p . 



The computation of an L 2 -defect bound 5\ satisfying 

|| - Auj x - \uj\ - |wa| 3 ||l2 < 5\ (6.3) 
amounts to the computation of an integral over Q. 

Due to [TTJ every solution of (12. ip is symmetric with respect to reflection at the axes 
x\ = \ and X2 = \- Therefore it is useful to look for approximate solutions of the form 

N 

cu\(xx,x 2 ) = ctij sin(i7rxi) sin(j7rx 2 ). (6.4) 

i,j odd 

Using sum formulas for sin and cos one obtains for all n G No, £ G R 

sin((2n + \)txx) = ^2 ^ cos(2A;7ra;) + 1^ sin 
and thus u>\ can be written as follows: 



TlX ) 



k=l 



uj\{xi 1 X2) = «n sin(7rxi) sin(7ra;2) + 

I iv-i I 

L— J / k , 

&2k+i,2l+i I 2 cos(2i7rxi) + 1 J I 2 cos(2j7rx 2 ) + 1 J sin(7rxi) sin(7rx 2 ). (6.5) 



k,l=l \ i=l / \ j=l 
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Since cos(x) ranges in [—1,1] and sin(7ra;i) sin(7rx 2 ) is positive for (x 1 ,x 2 ) G £1 = (0, l) 2 , 
u\ will be positive if 

«u+ ^ a 2fe+ i, 2i+ i([-2fc + l,2fc + l])([-2Z + l,2Z + l]) C (0,oo). (6.6) 

k,l=l 

Condition ( 16. 6 p can easily be checked using interval arithmetic and is indeed always sat- 
isfied for our approximate solutions, since an turns out to be "dominant" and the higher 
coefficients decay quickly. Hence u\ is positive and one can omit the modulus in the 
computations. Therefore the integral in (I6.3P can be computed in closed form, since only 
products of trigonometric functions occur in the integrand. After calculating them, var- 
ious sums Ylf=i remain to be evaluated. In order to obtain a rigorous bound 5\, these 
computations (in contrast to those for obtaining uj\ as described above) need to be carried 
out in interval arithmetic [T3J, EZj, to take rounding errors into account. 
Note that the complexity in the evaluation of the defect integral in (16. 3p . without any 
further modifications, is 0(N 12 ) due to the term uf. Using some trick, it is however 
possible to reduce the complexity to 0(iV 6 ): 

Applying the sum formulas sin(a) sin(6) = i[cos(a — b) — cos(a + b)] and cos(a) cos(6) = 
|[cos(a — b) + cos(a + b)] one obtains: 



sva.yii'KX) sm^ 2 7rxj smyi^KX) sm.{i^KX) smyigKx) sm(^ 6 7ra;j = 

' TTX) 



~ qo S a 2 a 3 a 4 a 5 a 6 cos ((h + a 2 i2 + & 3 i 3 + a 4 u + a 5 i 5 + a e i 6 ) 



32 

^2,^3,^4, 
0"5,0"6G{-1,1} 



Since £ cos(nvrx) dx = j J for ^ £ Z\{0} } =: 6n > We get 
J u x (x 1 ,x 2 ) (i d(x 1 ,x 2 ) = ^4 Yl & 2 - . . . ■ a e - p 2 - . . . ■ p 6 - 

n ff2,-^66{-l,l} p2,.--,P6e{-l,l} 

N N 

^ l ^ ] ^il+o-2i2 + ---+<J6i6^jl+p2j2 + ---+P6j6 a iljl ' • • • ' a i 

h,— 1*6=1 3l>-j6=l 

Setting aty := for G Z 2 \{1, . . . , N} 2 the previous sum can be rewritten as 



1024 



E 



a 2 ■ . . . ■ (7 6 • p 2 ■ . . . ■ p 6 - 



f>2 — .P66{-1,1} 



3N 



3N 



E E 

=-27V+l l=-2N+l 



( 

E E 

\+0-3«3=fc +P3J3=/ 



a hji a i2j2 a hji 



E E 



a j4j4 a *5j5 a *6j6 



. . (T4«4+o-5i5 P4.74+P5.75 
/ \+cr 6 i 6 = -k +p 6 j 6 = -Z 



For fixed o~i,Pi, k and Z each of the two double-sums in parentheses is 0(N 4 ). Since they 
are independent, the product is still 0(N 4 ). The sums over k and I then give 0(N 6 ), 
whereas the sums over <jj and pi do not change the complexity. 
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x is only 






v^3V 


if 


0~2 


l^k=2-N 


if 


<7 2 


l^k=-2N+l 


if 


Go 



Moreover the sum Ylt—- 

= 1,(73 = 1 
0-3 = -1 

= -1,(73 = -!. 

Similarly, also certain constellations of (74, (75, oe reduce the fc-sum, and of course analogous 
reductions are possible for the /-sum. Since ay = if z or j is even, the result does not 
change if the sum is only taken over odd values of i n ,j n , k and I. 

Remark 6.1. a) Computing trigonometric sums in an efficient way is an object of inves- 
tigation since a very long time, but up to our knowledge the above complexity reduction 
has not been published before. 

b) As an alternative to the closed form integration described above, we also tried quad- 
rature for computing the defect integral, but due to the necessity of computing a safe 
remainder term bound in this case, we ended up in a very high numerical effort, since a 
large number of quadrature points had to be chosen. So practically closed-form integra- 
tion turned out to be more efficient, although its complexity (as iV — > 00) is higher than 
the quadrature complexity. 

Once an L 2 -defect bound 5\ (satisfying (16.31) ) has been computed, an if _1 -defect bound 
5\ (satisfying (13.11) ) is easily obtained via the embedding 

< ^=L=NU 2 (ueL 2 (Q)) (6.7) 

which is a result of the corresponding dual embedding (15. 9p . Indeed, (I6.3P and (16. 7p imply 
that 

5\ := — , =5\ 

satisfies (13. ip . 

The estimate (16 .7p is suboptimal but, under practical aspects, seems to be the most 
suitable way for obtaining an if _1 -bound for the defect. At this point we also wish to 
remark that, as an alternative to the weak solutions approach used in this paper, we 
could also have aimed at a computer-assisted proof for strong solutions (see [23]), leading 
to if 2 -and C°-error bounds; in this case an L 2 -bound is needed directly (rather than an 
if ^-bound). 

7. Computation of K\ for fixed A 

For computing a constant K\ satisfying (13. 2p . we use the isometric isomorphism 

[ u ^ —Am + au J v ' 

and note that Q~ X L(\^ X ) : H\(^X) — > Hq(Q) is (•, • ^-symmetric since 

($ _1 L( A)a , A )[it],i'} i7 i = / [Vu ■ Vv - Xuv - 3\u x \u\uv} alx, (7.2) 

Jn 

and hence selfadjoint. Since ||£(a,u> a )M = ll ( ^ ) ~ 1 -^(A,w A )MII-f/i) (13. 2p thus holds for any 

K\ > [min : /i is in the spectrum of $ _1 L(a, Wa )}] 1 , (7.3) 
provided the min is positive. 
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A particular consequence of (17. 2p is that 

((J - S- 1 /^,) [u],u) B i = [ W x u 2 dx (u G H^Q)) (7.4) 

Jn 

where 

W x (x) := a + \ + 3\ux(x)\u x (x) (i6fl). (7.5) 

Note that, due to the positivity of our approximate solutions ux established in Section 
6, the modulus can be omitted here, which again facilitates numerical computations. 
Choosing a positive parameter a in the i^Q-product (recall that we actually chose a := 1), 
we obtain W\ > on Cl. Thus, (I7.4p shows that all eigenvalues \i of $ -1 L( AjWA ) are less 
than 1, and that its essential spectrum consists of the single point 1. Therefore, (17. 3p 
requires the computation of eigenvalue bounds for the eigenvalue (s) \i neighboring 0. 
Using the transformation k = 1/(1 — //), the eigenvalue problem $~ 1 L(a, Wa )[ , u] = \iu is 
easily seen to be equivalent to 



—Au + ou = kW\u, 



or, in weak formulation, 



(u,v) H i = k / Wxuvdx (v G Hq(Q)), (7.6) 
Jn 

and we are interested in bounds to the eigenvalue(s) k neighboring 1. It is therefore 
sufficient to compute two-sided bounds to the first m eigenvalues K\ < ■ • • < K m of 
problem f !7.6p . where m is (at least) such that K m > 1. In all our practical examples, 
the computed enclosures such that Ki < 1 < k 2 , whence by (17.31) and 

k = 1/(1 — //) we can choose 

^:=max{-^-,^-}. (7.7) 

Remark 7.1. By [11] and the fact that ux is symmetric with respect to reflection at 
the axes x\ — \ and X2 = |, all occurring function spaces can be replaced by their 
intersection with the class of reflection symmetric functions. This has the advantage that 
some eigenvalues Ki drop out, which possibly reduces the constant K x . 



The desired eigenvalue bounds for problem (17.61) can be obtained by computer-assisted 
means of their own. For example, upper bounds to Ki, . . . , K m (with m G N given) are 
easily and efficiently computed by the Rayleigh-Ritz method [26] : 

Let 0i, ... , (p m G Hq{VL) denote linearly independent trial functions, for example approx- 
imate eigenfunctions obtained by numerical means, and form the matrices 

Ax := {{<Pi,<Pj) H i)i,j=i,...,m, A := I J Wx&fijdx J 

\Q ' i,j=l,...,m 

Then, with Ai < • • • < A m denoting the eigenvalues of the matrix eigenvalue problem 

A\x = AAox 

(which can be enclosed by means of verifying numerical linear algebra; see [3]), the 
Rayleigh-Ritz method gives 

K{ < Aj for i — 1, . . . , m. 
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However, also lower eigenvalue bounds are needed, which constitute a more complicated 
task than upper bounds. The most accurate method for this purpose has been proposed 
by Lehmann [15], and improved by Goerisch concerning its range of applicability [I]. Its 
numerical core is again (as in the Rayleigh-Ritz method) a matrix eigenvalue problem, 
but the accompanying analysis is more involved. In particular, in order to compute lower 
bounds to the first m eigenvalues, a rough lower bound to the (m + l)-st eigenvalue must 
be known already. This a priori information can usually be obtained via a homotopy 
method connecting a simple "base problem" with known eigenvalues to the given eigen- 
value problem, such that all eigenvalues increase (index-wise) along the homotopy; see 
[22] or [5] for details on this method, a detailed description of which would be beyond the 
scope of this article. In fact, [5] contains the newest version of the homotopy method, 
where only very small (2 x 2 or even lxl) matrix eigenvalue problems need to be treated 
rigorously in the course of the homotopy. 

Finding a base problem for problem (I7.6p . and a suitable homotopy connecting them, 
is rather simple here since Q is a bounded rectangle, whence the eigenvalues of —A on 
Hq(Q) are known: We choose a constant upper bound c for |wa|wa = ^ x on Q, and the 
coefficient homotopy 

W ( x s \x) := a + A + 3[(1 - s)c + sco x (x) 2 ] (x E tt, < s < 1). 

Then, the family of eigenvalue problems 

-Am + an = k {s) W ( x s) u 

connects the explicitly solvable constant-coefficient base problem (s = 0) to problem (17. 6p 
(s = 1), and the eigenvalues increase in s, since the Rayleigh quotient does, by Poincare's 
min-max principle. 



8. Computation of branches (cu\), (S\), (K\), (a x ) 

In the previous section we described how to compute approximations u x for a grid of 
finitely many values of A within [0, Ai(fi)). After selecting A (among these) according to 
Section 5, we are left with a grid 

= A < A 1 < ■ ■ • < A A/ = A 

and approximate solutions u % = uj x * E Hq(D,) fl L°°(Q) (i = 0,...,M). Furthermore, 
according to the methods described in the previous sections, we can compute bounds 
6* = 8 X i and K* = K X i such that (EE]) and hold at A = A\ 

Now we define a piecewise linear (and hence continuous) approximate solution branch 
([0,~X]^H^n),X^co x ) by 



u x := ^^T^- 1 + A_ (A «-i < A < X, i = 1, . . . , M). (8.1) 



To compute corresponding defect bounds 5 X , we fix i E {1, . . . , M} and A E [A* 1 , A 1 ], and 
let t := (A - X i ~ 1 )/(X i - A^ 1 ) E [0, 1], whence 

A = (1 - t)^- 1 + tX\ u x = (1 - tp 1 - 1 + tu\ (8.2) 
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Using the classical linear interpolation error bound we obtain, for fixed iefl, 
\u x {x) 3 - [(1 - t)J~\xf + tJ(xf] I < 



< - max 

2 se[o,i] 



^[(l-s)oo i - 1 (x) + soo i (x 13 



t(l-t) 



< - max [(1 - s^-^x) + sJ(x)] ■ (u\x) - J- X {x)f 



4 se[o,i] 
3 

< -maxjUw^^loo, ||w l ||oo} \\u) 1 -u 



i-l||2 



13) 



\Xoo x (x) - [(1 - QX-W-^x) + tAV(x)] 



< - max 

2 se[o,i] 



d 



— [((1 - s)X i ~ 1 + SX)((1 - S)^-\X) + «,<(*))] 



t(l-t) 



< -(A 1 - A*- 1 )!^-^*- 11 



U) 



Since < Ci||w||oo for all w G L°°(f2), with Ci denoting an embedding constant for 

the embedding H$(Q) (e.g. Cj = v1fi[C 2 ), (ESD and (JH3D imply 



< -d maxjUa; 4 - 1 1 



oo, ll^lUjlla; 1 -J Yoc 



-■■Pi, (8.5) 



1 

—i 

4" 



|Aw A - [(1 - *)A t- V _1 + *AV]|| H -i < - A Ci(\ l - A* -1 )^ - w^Hqo =: n. (8.6) 



Now (JH2D, (JUSD, dHSD give 
II - Au\ - \uj\ - uIWh-i 



< (1 - t)\\ - Aw'" 1 - A'-V- 1 - (u/ -1 )iff-i + t|| - Aw 1 - AV - (u/)l^-i + r< + p t 

< maxitf- 1 , 5*} + n + p { =: 8 X . (8.7) 

Thus, we obtain a branch (<5a)a€[o,a] °f defect bounds which is constant on each sub interval 
[A* -1 , A*]. In the points A 1 , . . . , A M_1 , 8\ is possibly doubly defined by (18. 7p . in which case 
we choose the smaller of the two values. Hence, ([0,A] -4 E,A 4 5\) is lower semi- 
continuous. 

Note that 5\ given by (18 .7p is "small" if and 5 l are small (i.e. if the approximations 
u/ -1 and oj l have been computed with sufficient accuracy; see Remark 3.1a)) and if Pi,Ti 
are small (i.e. if the grid is chosen sufficiently fine; see (I8.5p . (I8.6P ). 



In order to compute bounds K\ satisfying (13. 2p for A G [0, A], with u\ given by ( 18. ip . we 
fix % G {1, . . . , M - 1} and A G + A*), §(A* + A m )] . Then, 

|A - A* | < - max{A l - A 4 " 1 , A m - X 1 } =: p h 

II^A-^llifi < -max{||u/ - w l ^^i, ||w l+1 - w l || H i} =: i/ f , 
whence a coefficient perturbation result given in pJH Lemma 3.2] implies: If 
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d ■= K l 



1 



Hi + 27{\\uj l \\ L A + C i v i )v i < 1, 



(8.9) 



Ai(fi) + a 



then (E2D holds for 



1-0' 



(8.10) 



Note that (18. 9p is indeed satisfied if the grid is chosen sufficiently fine, since then /i; and 



where K\ is possibly doubly defined by ( I8.10p . we obtain a lower semi-continuous, piece- 
wise constant branch ([0, A] — > R, A i— >■ K\). 

According to the above construction, both A H- 5 A and A i-> K\ are constant on the 2M 
half-intervals. Moreover, ( 18.1 j) implies that, for % = 1, . . . , M, 



and again we choose the smaller of the two values at the points of double definition. 
Using these bounds, the crucial inequalities ( 13 .4p and (13. 5p (which have to be satisfied for 
all A e [0, A]) result in finitely many inequalities which can be fulfilled with "small" and 
piecewise constant a>\ if 5\ is sufficiently small, i.e. if u°, . . . ,u M have been computed 
with sufficient accuracy (see Remark 3.1a)) and if the grid has been chosen sufficiently 
fine (see (18.51) - (18. 7p ). Moreover, since A 5\, A K\ and the above piecewise constant 
upper bound for ||wa||l 4 are lower semi-continuous, the structure of the inequalities (13 .4p 
and (13. 5p clearly shows that also A a\ can be chosen to be lower semi-continuous, 
as required in Section 4. Finally, since (I3.5P now consists in fact of finitely many strict 
inequalities, a uniform (A- independent) rj > can be chosen in Theorem 3.1b), as needed 
for Theorem 4.1. 



All computations have been performed on an AMD Athlon Dual Core 4800+ (2.4GHz) 
processor, using MATLAB (version R2010a) and the interval toolbox INTLAB [27]. For 
some of the time consuming nested sums occurring in the computations, we used moreover 
mexfunctions to outsource these calculations to C++. For these parts of the program we 
used C-XSC [13] to verify the results. Our source code can be found on our webpage^. 
In the following, we report on some more detailed numerical results. 



http : //www. math. kit . edu/iana2/~roth/page/publ/en 



Ut are "small" by <KE\f . 




wa||l-i < 



maxdlw 4 - 1 !!^, Kllw*" 1 !!^ + |M| L 4)} for A e [A*" 1 , ^A^ 1 + A*)] 
max{|(||a; < - 1 |U* + \\u>%*), ||w*M} for A E [|(A i_1 + A J ), A*] 



9. Numerical results 
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Figure 1. Curve (A, ||wa||c») with samples of u\ in the case p = 3 



Using A = 18.5 (which is not the minimally possible choice; e.g. A = 15.7 could have 
been chosen) and M + 1 = 94 values = A < A 1 < ■ • • < A 93 = 18.5 (with A 1 = 
0.1, A 2 = 0.3 and the remaining gridpoints equally spaced with distance 0.2) we computed 
approximations u°, . . . , u 93 with N = 16 in (16.21) . as well as defect bounds 5°, . . . , 5 93 and 
constants K°, . . . , K 93 , by the methods described in Section 6 and 7. 
Figure 1 shows an approximate branch [0,27r 2 ) — >• R, A (->■ ||wa||oo- The continuous plot 
has been created by interpolation of the above grid points A- 7 , plus some more grid points 
between 18.5 and 27r 2 , where we computed additional approximations. 
For some selected values of A, Table 1 shows, with an obvious sub- and superscript no- 
tation for enclosing intervals, the computed eigenvalue bounds for problem (17. 6p (giving 
K\ by (17. 7p ). These were obtained using the Rayleigh-Ritz and the Lehmann-Goerisch 
method, and the homotopy method briefly mentioned at the end of Section 7 (exploiting 
also the symmetry considerations addressed in Remark 7.1). The integer m, needed for 
these procedures, has been chosen different (between 3 and 10) for different values of A, 
according to the outcome of the homotopy. This resulted in a slightly different quality of 
the eigenvalue enclosures. 
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«2 




0.34350814513^ 


2.492570^ 


LU2.7 


0.37521912233l|g 


2.6221837^1 


w 6.7 


0.4373273950^ 


2.87378161^ 


^10.7 


0.52752354636^1 


3.223417042^ 


^14.7 


0.6676848259^ 


3.725209290H 


^18.5 


0.89237445994^1 


4.46288110^1 



Table 1. Eigenvalue enclosures for the first two eigenvalues 

Table 2 contains, for some selected of the 186 A-half-intervals, 

a) the defect bounds 5\ obtained by (18. 7p from the grid-point defect bounds <5 l_1 ,<5 1 , 
and from the grid- width characteristics pi,Ti defined in (18. 5p . (18. 6p . 

b) the constants K\ obtained by (18.101) from the grid-point constants K l and the 
grid- width parameters defined in (18. 8p (note that /ij = 0.1 for all i), 

c) the error bounds a\ computed according to (13 .4p . (13. 5p . 

Thus, Corolllary I2.1[ together with all the considerations in the previous sections, proves 
Theorem 11.11 



A-interval 


h 


K x 


a x 


[0,0.05) 


0.0005943 


1.7443526 


0.0010378 


(2,2.1) 


0.0023344 


1.7707941 


0.0041521 


(6,6.1) 


0.0022937 


1.6669879 


0.0038369 


(10,10.1) 


0.0023644 


1.5677657 


0.0037168 


(14,14.1) 


0.0026980 


1.9582604 


0.0053028 


(16,16.1) 


0.0031531 


3.2267762 


0.0102701 


(18.4,18.5] 


0.0050056 


13.8930543 


0.0882899 



Table 2. 

Acknowledgments. The authors are grateful to two anonymous referees for their helpful 
remarks and suggestions. 
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